* Mechanisms analysis: cap format X accessibility of food stores
* 3/9/2025

clear all
set more off

program main
	state_fips_xwalk
	read_usda_data
	merge_tract_puma
	prepare_acs
	
	geo_regressions, cap_format(AE) low_access_var(mean_one) // low access at 1 and 10 mi.
	geo_regressions, cap_format(standard) low_access_var(mean_one)
end

program state_fips_xwalk
	use "$for_analysis/state_for_regressions_actual_ssi", clear
	keep statefip
	duplicates drop
	
	decode statefip, gen(State)
	replace State = strproper(State)
	save "$intermediate/state_fips_xwalk", replace
end

program read_usda_data // clean and save USDA food access dataset
	syntax[, cap_only(str)]
	
	// downloaded from https://www.ers.usda.gov/data-products/food-access-research-atlas/download-the-data
	// most recently downloaded May 26, 2025
	import excel using "$raw/USDA_food_access_2015/FoodAccessResearchAtlasData2015.xlsx", ///
		sheet("Food Access Research Atlas") clear firstrow
	drop Tract* LILA*

	if "`cap_only'" == "true" { // restrict to treated states
		keep if inlist(State, "South Carolina", "Mississippi", "Washington", "Texas", "New York") | ///
			inlist(State, "Massachusetts", "Florida", "North Carolina", "Pennsylvania", "Virginia") | ///
			inlist(State, "Kentucky", "Louisiana", "Michigan", "New Jersey", "Arizona", "South Dakota") | ///
			inlist(State, "New Mexico", "Maryland")
	}
	
	replace State = "District Of Columbia" if State == "District of Columbia"
	merge m:1 State using "$intermediate/state_fips_xwalk", assert(3) keep(3) nogen
	
	gisid CensusTract
	save "$intermediate/usda_access", replace
end

program merge_tract_puma  // merge in xwalk from 2010 census tract to 2000 PUMA
	* produced from "Geocorr 2018: Geographic Correspondence Engine", at https://mcdc.missouri.edu/applications/geocorr2018.html
	* converts 2010 census tracts to (1) 2000 PUMAs and (2) 2012 PUMAs; I read this file as a .rtf and converted to .xls in Excel
	* note: when read into Excel, make sure all variables are stored as "Text", not "General"
	import excel "$raw/geo_data/puma_tract_xwalk.xlsx", clear firstrow

	replace tract = subinstr(tract, ".", "", .) 
	tostring county, replace
	gen CensusTract = county + tract
	ren state statefip
	destring statefip pop10, replace
	drop *name

	merge m:1 statefip CensusTract using "$intermediate/usda_access", assert(1 2 3) gen(merge_usda_xwalk)
	tab State if merge_usda_xwalk == 2  // a few South Dakota and Alaska census tracts did not merge
	tab statefip if merge_usda_xwalk == 2
	keep if merge_usda_xwalk == 3
	cap drop merge_usda_xwalk

	* calculate average measure of food access by PUMA
	keep statefip puma2k puma12 County pop10 CensusTract LA1and10 ///
		LAhalfand10 LA1and20 LATractsVehicle_20 laseniorshalfshare
	
	foreach geo of varlist puma2k puma12 {  // both 2000 and 2012 PUMA
		preserve
			* weight by tract population
			bys statefip `geo': egen total_pop = sum(pop10)
			gen pop_weight = pop10/total_pop
			
			gcollapse (mean) LA1and10 LAhalfand10 LA1and20 LATractsVehicle_20 ///
				laseniorshalfshare [aw = pop_weight], by(statefip `geo')

			ren `geo' puma
			if "`geo'" == "puma2k" {
				gen merge_year = 2000
			}
			else {
				gen merge_year = 2012
			}
			save "$intermediate/`geo'_food_access", replace
		restore
	}
	
	* create dataset for map in R, using 2000 PUMAs
	use "$intermediate/puma2k_food_access", clear 
	keep if inlist(statefip, 45, 12, 25, 36, 48, 53, 28, 37) | ///
		inlist(statefip, 42, 51, 21, 22, 26, 46, 4, 34, 35, 24) // restrict to CAP states
	tostring statefip, replace
	tostring puma, replace
	replace statefip = "0" + statefip if strlen(statefip) < 2
	assert strlen(puma) == 5
	assert strlen(statefip) == 2
	order statefip puma
	ren (statefip puma) (STATEFIP PUMA)
	save "$intermediate/puma_food_access_R", replace
	export delimited "$intermediate/puma_food_access.csv", replace
end

program prepare_acs
	use "$for_analysis/state_for_regressions_actual_ssi", clear
	drop if year < 2005  // geographic data recorded in 2005 and later
	assert !mi(puma)
	gen merge_year = 2000 if year < 2012
	replace merge_year = 2012 if year >= 2012 // ACS data uses 2010 PUMAs for 2012 and later
	tostring puma, replace
	
	replace puma = "0" + puma if strlen(puma) == 4
	replace puma = "00" + puma if strlen(puma) == 3
	assert strlen(puma) == 5
	
	merge m:1 statefip puma merge_year using "$intermediate/puma2k_food_access", ///
		assert(1 2 3) keep(1 3) gen(merge2000) // 2000 PUMAs, for years 2005-2011 of ACS
	merge m:1 statefip puma merge_year using "$intermediate/puma12_food_access", ///
		assert(1 2 3 4) keep(1 3 4) update gen(merge2012) // 2010 PUMAs, for years 2012 on of ACS
		
	tab statefip if merge2012 == 1 & merge2000 == 1 // just one PUMA in Louisiana is problematic
	keep if merge2000 == 3 | merge2012 == 3
	cap drop merge*
	
	ren (LA1and10 LAhalfand10 LATractsVehicle_20 laseniorshalfshare) ///
		(mean_one mean_half mean_vcle mean_senior)
	assert !mi(mean_one)
	save "$for_analysis/state_geo_for_regressions", replace
end

program geo_regressions
	syntax, cap_format(str) low_access_var(str)

	use "$for_analysis/state_geo_for_regressions", clear
	if "`cap_format'" == "standard" {
		keep if standard == 1 | mi(startyear)
	}
	else if "`cap_format'" == "AE" {
		keep if inlist(statefip, 36, 25, 42) | mi(startyear) // AE states plus all controls
	}
	else if "`cap_format'" == "modified" {
		keep if standard == 0 | mi(startyear)
	}
	
	local indiv_controls = "i.sex i.race i.hispan i.edcat age i.diffphys"
	local state_snap_controls = "i.has_bbce i.has_call i.has_faceini i.has_facerec i.has_oapp i.has_no_fp"

	su `low_access_var', d
	gen low_access = (`low_access_var' >= `r(p50)')

	* no controls
	reghdfe snap i.inter##i.low_access [w = perwt], ///
		a(i.statefip i.year) vce(cl statefip)
	est sto reg1
	qui estadd local cont " ", replace
	qui estadd local unemp " ", replace
	qui estadd local hhtype "\checkmark", replace
	qui estadd local yearFE "\checkmark", replace
	
	* individual controls
	reghdfe snap i.inter##i.low_access `indiv_controls' [w = perwt], ///
		a(i.statefip i.year) vce(cl statefip)
	est sto reg2
	qui estadd local cont "\checkmark", replace
	qui estadd local unemp " ", replace
	qui estadd local hhtype "\checkmark", replace
	qui estadd local yearFE "\checkmark", replace
	
	* individual & state controls
	reghdfe snap i.inter##i.low_access ///
		`indiv_controls' state_unemp `state_snap_controls' [w = perwt], ///
		a(i.statefip i.year) vce(cl statefip)
	est sto reg3
	qui estadd local cont "\checkmark", replace
	qui estadd local unemp "\checkmark", replace
	qui estadd local hhtype "\checkmark", replace
	qui estadd local yearFE "\checkmark", replace
	
	* generate table
	esttab reg1 reg2 reg3, ///
		keep(1.inter 1.inter#1.low_access) order(1.inter 1.inter#1.low_access) ///
		varlabels(1.inter "Treated x Post" ///
		1.inter#1.low_access "Treated x Post x [Low Access to SNAP Retailers]") ///
		se(%5.3f) star(* 0.10 ** 0.05 *** 0.01) stats(yearFE hhtype cont unemp N, ///
		labels("Year FE" "State FE" "Individual Controls" "State Controls" "N" ) ///
		fmt("%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc")) scalars() ///
		nomtitles label nonote booktabs
	
	esttab reg1 reg2 reg3 ///
		using "$output/twfe_interact_low_access_`cap_format'", replace ///
		keep(1.inter 1.inter#1.low_access) order(1.inter 1.inter#1.low_access) ///
		varlabels(1.inter "Treated x Post" ///
		1.inter#1.low_access "Treated x Post x [Low Access to SNAP Retailers]") ///
		se(%5.3f) star(* 0.10 ** 0.05 *** 0.01) stats(yearFE hhtype cont unemp N, ///
		labels("Year FE" "State FE" "Individual Controls" "State Controls" "N" ) ///
		fmt("%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc")) scalars() ///
		nomtitles label nonote booktabs
end

* Execute
main
